Cell-specific models of hiPSC-CMs developed by the gradient-based parameter optimization method fitting two different action potential waveforms

Parameter optimization (PO) methods to determine the ionic current composition of experimental cardiac action potential (AP) waveform have been developed using a computer model of cardiac membrane excitation. However, it was suggested that fitting a single AP record in the PO method was not always successful in providing a unique answer because of a shortage of information. We found that the PO method worked perfectly if the PO method was applied to a pair of a control AP and a model output AP in which a single ionic current out of six current species, such as IKr, ICaL, INa, IKs, IKur or IbNSC was partially blocked in silico. When the target was replaced by a pair of experimental control and IKr-blocked records of APs generated spontaneously in human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs), the simultaneous fitting of the two waveforms by the PO method was hampered to some extent by the irregular slow fluctuations in the Vm recording and/or sporadic alteration in AP configurations in the hiPSC-CMs. This technical problem was largely removed by selecting stable segments of the records for the PO method. Moreover, the PO method was made fail-proof by running iteratively in identifying the optimized parameter set to reconstruct both the control and the IKr-blocked AP waveforms. In the lead potential analysis, the quantitative ionic mechanisms deduced from the optimized parameter set were totally consistent with the qualitative view of ionic mechanisms of AP so far described in physiological literature.


Cell-specific models of hiPSC-CMs developed by the gradient-based parameter optimization method fitting two different action potential waveforms
Yixin Zhang 1 , Futoshi Toyoda 2,3 , Yukiko Himeno 1,2* , Akinori Noma 1 & Akira Amano 1 Parameter optimization (PO) methods to determine the ionic current composition of experimental cardiac action potential (AP) waveform have been developed using a computer model of cardiac membrane excitation.However, it was suggested that fitting a single AP record in the PO method was not always successful in providing a unique answer because of a shortage of information.We found that the PO method worked perfectly if the PO method was applied to a pair of a control AP and a model output AP in which a single ionic current out of six current species, such as I Kr , I CaL , I Na , I Ks , I Kur or I bNSC was partially blocked in silico.When the target was replaced by a pair of experimental control and I Kr -blocked records of APs generated spontaneously in human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs), the simultaneous fitting of the two waveforms by the PO method was hampered to some extent by the irregular slow fluctuations in the V m recording and/or sporadic alteration in AP configurations in the hiPSC-CMs.This technical problem was largely removed by selecting stable segments of the records for the PO method.Moreover, the PO method was made failproof by running iteratively in identifying the optimized parameter set to reconstruct both the control and the I Kr -blocked AP waveforms.In the lead potential analysis, the quantitative ionic mechanisms deduced from the optimized parameter set were totally consistent with the qualitative view of ionic mechanisms of AP so far described in physiological literature.
The AP signal was filtered at 5 kHz and digitized at 20 kHz.A selective I Kr blocker, E-4031 (Wako, Japan) was dissolved in distilled water to yield 5 mM stock solution and added to the bath solution at the concentrations of 0.01-1 μM.

The baseline model of hiPSC-CM membrane excitation
The baseline model of hiPSC-CMs is a variant of the human ventricular cell model (hVC model) adapted to experimental measurements in hiPSC-CMs 25 .The hVC model has been fully described in literature 9,10 and shares many comparable characteristics with other human models so far published 7,8 .The model structure of the hVC model consists of the cell membrane with fifteen ionic channel species and two ion exchangers, the sarcoplasmic reticulum equipped with the refined Ca 2+ releasing units coupled with the L-type Ca 2+ channels on the cell membrane at the microscopic dyadic space, the contractile fibers, and the cytosolic Ca 2+ diffusion spaces containing several Ca 2+ -binding proteins.
The source code of the present hiPSC-CM model was written in Visual Basic and is available from our e-Heart website (http:// www.ehear tsim.com/ en/ downl oads/), and the model equations are described in the supplemental materials in the present study.
In the present study, the net membrane ion current (I tot_cell ) is calculated as the sum of nine ion channel currents and two ion transporters (I NaK and I NCX ) (Eq. 1).
A few ion channels listed in Eq. (1) or in supplemental materials were not necessarily expressed in every cell to a significant level to affect the AP waveform.In such cases, we removed the current from the initial set of PO parameters through the initial manual fitting process explained in the last part of this section.
The I Na is composed of two components; one is the transient type (I NaT ) and the other shows the very slow inactivation (I NaL ).I NaL plays a significant role in both slow diastolic depolarization as well as in the AP plateau.In the present study, the Na + mediated sustained inward current, I st [26][27][28] is included in the I CaL for simplicity.Rationale for this simplification was derived from the references suggesting Cav1.3 as a molecular determinant of I st 29,30 .We excluded background currents of much smaller amplitude, such as I KACh , I KATP , I (lCa) , and I Cab from the parameter optimization and adjusted only the background non-selective cation current (I bNSC ) of significant amplitude for the sake of simplicity [31][32][33] .The I bNSC was re-defined in the present study as a time-independent net current, which remained after blocking all time-dependent currents.
The parameter optimization in the present study was conducted under the assumption that the variable AP configurations in the hiPSC-CMs were caused by different expression levels of ion channels on the cell surface membrane, while the channel open-and close-state transition kinetics remain the same as described in the hVC cell model.Before starting the PO method, we manually adjusted the model parameters (current magnitude of individual currents in Eq. 1) to fit the model output to the target AP waveform on the graphic display.Usually, the 9 ionic channel currents were all examined for their capability to simulate the experimental AP configuration by the model output, and those without an obvious improvement of the model output were excluded from the PO parameters.This selection of the optimizing parameters contributes to reducing the computational time of the PO method.The computational time increases exponentially when the number of the optimizing parameters increase and also if the parameter which does not affect the AP waveform was included.Note that, the number of the optimizing parameters to perform the PO method in the realistic computational time was around six with the computational resource used in this study.

A brief note on membrane excitation
The membrane excitation of the model is generated through charging and discharging the membrane capacitance (C m ) by the net ionic current (I tot_cell ) across the cell membrane (Eq.2).The driving force for the ionic current is provided by the potential difference between V m and the equilibrium potential (E x ) (Eq. 3).The net electrical conductance of the channel is changed by the dynamic changes in the open probability (pO) of the channel, which is mostly V m -dependent through the V m -dependent rate constants ( α, β ) of the opening and closing conforma- tion changes of the channel (Eqs.4, 5).
(1) I tot_cell = I Na + I CaL + I ha + I K1 + I Kr + I Ks + I Kur + I Kto + I bNSC + I NaK + I NCX (2) For the parameter optimization, the limiting conductance G x of each ionic current (Eq.3) is represented by a product of I xmag and I xsf .The former is a reference magnitude of conductance and the latter is the scaling factor (0 < I xsf ).
When the constant field (GHK) equation is used, the current density amplitude I x is given by, where P x is the limiting permeability of ion x across the cell membrane and S i and S o are the internal and external ion concentrations, respectively.The unit of P x is pA/pF/mM across the membrane.The I xsf is also used for the PO parameter.

Pattern search method; MM test, EM test, and multi-run test
For a system showing a relatively simple gradient of MSE along the parameter axis, gradient-based optimization methods are more efficient in general than stochastic methods.We used one of the basic gradient-based optimization methods, the pattern search (PS) method [34][35][36] .The computer program code of the pattern search 37 is simple and does not require derivatives of the objective function.We implemented the program code into a homemade program for data analysis.
The parameter optimization is driven by calculating the error function which uses the mean square error (MSE) between the target AP waveform (V m,t ) and the output of the baseline model (V m,a ) (Eq. 11), where V m,a represents adaptive V m (the model output) generated by adjusting the randomized set of I xsf s in the baseline model.The target V m,t is the AP waveform of the intact baseline model.The magnitude of MSE was stabilized by obtaining a stable rhythm of spontaneous APs through repetitive numerical integration of the model, usually, 30 ~ 100 spontaneous cycles were repeated for a new set of I xsf s.The MSE was calculated within a time window which included a single cycle of the spontaneous AP generation.N is the number of digitized V m points with a time interval of 0.1 ms.
The present study used the weighted sum of MSE.The single cycle of the spontaneous AP was separated into three segments.The first segment starts from − 20 mV during the repolarizing phase to maximum diastolic potential (MDP) and the weight is 5 ~ 6.The second segment covers MDP, SDD, and the foot of AP discharge and the weight is 1, and the third segment covers the AP plateau until to repolarizing phase to − 20 mV with weight = 0.3 or 0.5.In the first segment, the magnitude of I Kr dominates to allow its determination by the PO method.The second segment is influenced mostly of I bNSC , I ha in addition to I NaL and I CaL , which takes a major role in generating the exponential depolarization in the late half of SDD.The third segment consists of the AP plateau, generated by I CaL , I NCX , I NaL , I Ks , and I Kur .By providing large weight to the first time window, I Krsf rapidly gets close to the optimal value independent of other parameters which means that the succeeding optimization process behaves as if the number of the optimizing parameters is one less.According to our experience, the preliminary manual fitting can be carried out smoothly in a sequence of these three segments.
All hiPSC-CMs used in the present study showed spontaneous discharge of AP.Therefore, to calculate MSE, it was essential to start all records from a particular V m value; − 20 mV was used in this study.We calculated the model until digitized V m crossed a range of − 20.000 mV > V m > − 20.004 mV.Thereafter, a record of a single AP cycle of the target or model output was obtained.For the graphic display, we usually connected two cycles of AP for better visualization.
In the PS method, a base point (BP) is determined by the initial parameter set.On both sides of BP, two new points (NP) are set by a ± step to calculate MSE for each NP, and the numerical integration of the hiPSC-CM model was conducted to obtain new AP waveforms.If MSE is improved, the NP is moved by the single step in the direction giving the smaller MSE.This pattern search is repeated for all dimensions of PO parameters one by one until the decrease of MSE saturates.Then all BPs are moved to the last NPs for all PO parameters, and simultaneously the step size is decreased.The movements of BP are repeated until the step size becomes smaller than a critical step size, which was preset by the user.( 4) The intracellular concentrations of Na + , Ca 2+ , and K + will change every time when the parameters are varied during the PS method.To minimize the variation of the ion concentration, the activities of both I NaK and I NCX were regulated.For I NaK , Eq. ( 12) was used.
where any deviation of [Na + ] (Na i ) from a standard [Na + ] (stdNa i ) was corrected by modifying crf NaK used for the negative feedback regulation in reference to a constant (orgcrf NaK ).The fbGain NaK was usually set at 0.004 for a stable regulation.The stdNa i = 5.8 mM and orgcrf NaK of 0.2606 were used.
The same format of equation was used to regulate the I NCX activity to keep the total of ionized and bound form of Ca (Ca content) within the cell, which is a sum of Ca content in the three Ca compartments and the sarcoplasmic reticulum (SR) ([Ca] t ).The fbGain NCX was 0.008.The std totCai = 68.5 attomole and orgcrf NCX = 2.4732.

Parameter randomization and the multi-run test
The final level of MSE obtained by a single cycle of the PS method fluctuates between different cycles.Therefore, the PS cycle was repeated for 200 ~ 300 cycles (multi-run test).The MSE as well as the model parameters were given as an average of the top 20 results of smaller MSE as described in Kohjitani et al. 25 .Every run of the PS method started with a randomized set of the PO parameters obtained by Eq. ( 14).R sf is given by a random function and is distributed between 0 ≤ R sf < 1 .The af is used to vary the ampli- tude of I xsf ; an af of 0.3 is usually used to randomize I xsf between 0.85 and 1.15.The time course of optimization was monitored by plotting I xsf in a two-dimensional parameter space on the graphic display; log(MSE) on the ordinate against I xsf on the abscissa.

The simultaneous fitting of two AP waveforms in the PS method
In the present PO method, the experimental data obtained under the specific block of I Kr using E-4031 in addition to the control AP waveform were used (see Doss et al. 38 for the I Kr block in hiPSC-CMs).Namely, MSE was calculated for both the control and the I Kr -blocked AP waveform with the model output, and the average of these two MSEs was used as the error function to drive the PO method.
The mechanisms underlying the modification of the cardiac AP waveform are implemented in the baseline model.Thus, the maximum diastolic potential, (MDP) as well as the slow diastolic depolarization (SDD) were shifted in the positive direction in most of hiPSC-CMs if G K (Eq. 3) was decreased by E-4031.Through these changes in V m , the V m -dependent open-close gatings were altered according to the gating kinetics to modify the AP waveform.Furthermore, under the specific I Kr block, a new constraint is applied to the PO; i.e. the G x of other I x should remain at the control magnitude.Thus, a common set of I xsf s except the I Krsf was used for both the control and I Kr -blocked AP waveforms.Note that, this constraint was also applied in the initial manual fitting process, thus the magnitude of I Kr block was manually estimated as the result in this process.

The lead potential (V L ) analysis
The final AP waveforms obtained by the PS method were examined for their validity by calculating the contributions of individual current components (I x s) to the AP waveform by using the lead potential (V L ) analysis 9,39 .Contributions of individual I x s are expressed in a unit of mV instead of pA/pF, where G x is the conductance (nS/pF) given by Eq. ( 3).The V m closely followed the time course of V L (Eq. 16) with the time constant τ (Eq.17).G x is the chord conductance in the case of currents described according to Ohm's law.In the case of constant field equations, the G x is obtained by differentiating GHK with respect to V m (Eqs.18 ~ 20) as the slope conductance (Eq.19), ( 12) The reversal potential (E x ) at a given V m was determined by Eq. ( 21).

Accuracy of the simultaneous two AP waveform fitting examined by the MM mode of the PO method
The graphic display of the MM multi-run tests conducted for each case of selective block of I Kr (60%), I Na (50%), I CaL (50%), I Ks (80%), and I bNSC (50%) were demonstrated in Fig. 1, where the movements of BPs in each run of PS method were plotted in the I xsf vs log(MSE) coordinates for each I xsf indicated at the top.The cell-specific model developed by the single AP waveform fitting was used as the baseline model.Although the PS parameters of the initial parameter set were randomized for each run of the PS method, the BPs were distributed within an obvious envelope and well converged at a normalized I xsf = 1 value in the baseline model.The minimum level of log(MSE) was mostly smaller than -2 (MSE < 0.01 mV 2 ), which is in the same range as those obtained in single AP waveform fitting (see Fig. 4 and Table 2 in Kohjitani et al. 25 ).In Fig. 2, the target AP waveforms (red) were superimposed by APs of the model output (black traces), which were the outcomes of the MM multi-run test of the PO method in Fig. 1.It is evident that the PO was nearly perfectly conducted when the log(MSE) below − 2 were obtained.All AP waveforms well represent the characteristic effect of the I x block.Namely, the blockage of the outward currents, I Kr and I Ks prolonged the APD, while the blockage of the inward currents, I CaL shortened the APD.In the case of I Na block, the SDD was markedly prolonged through the partial block of the I NaL component.During the prolonged SDD, the inactivation of I NaL might be partially removed to recover the APD.A similar prolongation of the APD was also observed by the I bNSC blockage.These findings are totally consistent with those expected from the physiological mechanisms of the AP generation.
Table 1 lists the mean and standard deviation of the magnitude of final I xsf s in the top 20 runs of PO in the MM tests performed to generate AP waveforms demonstrated in Fig. 2. The magnitudes of I xsf were normalized in reference to the standard values of the baseline model.In the case of blocked currents, the I xsf was normalized to the magnitude of blockage defined in each PS method.It is evident that the error of the I xsf is less than ~ 2% of the control and the SD of the mean was negligibly small, indicating that selective blockage of a given ionic current was correctly resolved by the PS method in the present study.
The level of [Na + ] i and [Ca] t indicate the final concentrations determined in the AP waveform under the selective I x blockage.The standard deviation (SD) of the mean was also negligibly small to support that the activity of the Na/K pump current and the Na/Ca exchange currents were well controlled to give a stable level of [Na + ] i and [Ca] t (Eqs.12 and 13) in the top 20 runs of the PS method.

Feasibility of applying the simultaneous two AP waveform fitting to experimental recordings of the selective I Kr block
According to the MM mode fitting, it is concluded that the algorithm of the PS method is relevant to get the cell-specific model using the two-waveform parameter fitting.However, it is asked if the EM mode is also feasible in analyzing the I xmag profile of experimental AP waveforms, which include various extra noise, such as the slow fluctuations of V m in the hiPSC-CMs.Furthermore, the EM fitting might be interfered with if the computer model equations failed to describe the intact mechanisms of the intact cardiac myocytes.This latter problem should be finally solved by performing a new experimental protocol, but the study of the EM mode PO method might suggest inappropriate equations of the model.Here, we examine the feasibility of applying the EM mode of the PS method to the experimental AP recordings.

Time course of the spontaneous APs after the application of E-4031 and results of the conventional PO method of using a single AP waveform
Three representative AP recordings in the hiPSC-CMs after exposure to an I Kr blocker E-4031 are demonstrated in Fig. 3A-C.The continuous recordings of V m in the hiPSC-CMs showed sporadic slow fluctuations of V m and/ or irregular cycle length of the spontaneous rhythm.In the record of a single myocyte (TC11) in Fig. 3A, for example, irregular fluctuations in MDP in addition to sporadic changes in AP configurations were obvious.The myocyte TC13 in Fig. 3B shows a relatively smooth time course of the MDP and OS change before and after the I Kr block.In the record of myocyte TC12 in Fig. 3C, the fluctuations in AP amplitude after about 12,000 ms were of unknown origin.These fluctuations were observed in most of the hiPSC-CMs.
Before conducting the new two-waveform fitting, we firstly examined the conventional PO method to apply the baseline model to individual AP waveforms.The baseline model was obtained by applying the PS method to an AP waveform during the control condition in the same myocytes.The E-4031 was applied by switching the perfusate using a cock, which connected the outlet tube to the reservoir of test solutions about 40 cm apart from the inlet of the recording chamber to avoid mechanical vibration.It was difficult to precisely indicate the onset time of the I Kr block by the E-4031 solution.However, the progressive positive shift of MDP induced by E-4031 was obvious.Another complication might be caused by the progress of the I Kr block during a single AP cycle (t 2 − t 1 ).Assuming an exponential development of the I Kr block, the progress was calculated using Eq. ( 22), in which br t2 indicates the block ratio at time t 2 .If a rate constant τ is assumed to be 30 s, the difference of br might be ~ 1% during the AP interval of 0.3 s. we assumed that an error due to this transitional effect might be negligibly small in testing the application of the PO method.A representative result of I xsf s optimized by the single waveform fitting of the PS method applied sequentially was illustrated for individual AP waveforms in TC12 in Fig. 3D.The I Krsf (light blue) was decreased smoothly in an exponential manner.The fluctuations in I Nasf (steel blue), I Kursf (yellow) and I bNSCsf (blue) were less than ~ 20%.However, I CaLsf (orange), largely deviated.This deviation of the I CaLsf might be attributed not only to the extra fluctuations in the continuous recording of the spontaneous activity, but also to the shortage of information provided by the experimental single AP waveform.

Results of applying the simultaneous two-AP-waveform fitting of the EM mode of the PS method.
From the limited length of the continuous AP record, stable AP waveforms of control and I Kr -blockage were selected for the EM mode of the simultaneous two-AP waveform fitting.The red traces in Fig. 3A-C indicated individual AP waveforms; c the reference and b1 and b2 the test waveforms used for the simultaneous fitting of pairs; c and b1, and c and b2 waveforms.Results of the PS method are shown in Fig. 4. Differently from the MM test of the PS method (Fig. 2), a slight gap remained after the end of the PO run between the target (red trace) and the model output (black trace), most probably due to the irregular slow fluctuations.Indeed, the final magnitude of log(MSE) is comparable to results obtained in the single AP wave fitting in EM mode in our previous study 25 .The result of the multi-run test shown in the right half of Fig. 4 showed clear convergence of the I Krsf , I CaLsf, and I bNSCsf , which gave the averages of top 20 results of smaller MSEs (vertical blue lines) mostly within the ± 5% Selective I Kr block resolved by the simultaneous fitting of two waveforms of the PS method.
The selective I x block has been proved in the MM mode of the PO method (Figs. 1, 2, and Table 1).In Table 1, the I xsf were all very close to 1, supporting the fine resolution of the selective I x block by the PS method.In the  EM mode of the PO method, however, the selective blocker which had a high selectivity to a specific channel was only available for I Kr .To prove the selective block of I Kr obtained in the EM mode PO, the magnitudes of I xsf s other than the I Krsf were compared in Fig. 5 at two different degrees of the I Kr block.The magnitudes of I Krsf were quite close to 1 in the multi-run test because the I Krmag was adjusted to the degree of the I Kr block, which was determined by manual fitting beforehand.For better visibility, the I xsf obtained with the higher I Kr block (green bars) was normalized in reference to that obtained at the lower I Kr block (blue bars in Fig. 5).It is evident that the deviation of I xsf was less than 10% in 16 measurements of I Krsf , and 10-20% in other 5 measurements.The large deviation of I ha in TC11 was partly attributed to the depolarized MDP level at the higher I Kr block, which will be discussed further in the next section.Considering the relatively large fluctuations of the experimental V m recordings (Fig. 3), it might be concluded that the EM mode of the PS method successfully revealed the selective I Kr block.

Correlation between two ion currents revealed by the EM test of PS method
When the error function MSE of two parameter sets are decreased to a minimum level in the PS method, the model output is almost perfectly fitted to the target AP waveform, which indicates that the I tot_cell in Eqs.(1) and ( 2) of two parameter sets are almost the same.If the two parameters I X1sf and I X2sf have ambiguity, these parameters in the final results of the PS method have multiple value combinations, however, they have almost From the above consideration, we can expect that if the number of the correlation of grade (a) is small, then the problem has small ambiguities.In Fig. 6, the grade of R 2 was represented by the thickness of the line, connecting each pair of I xsf s.It is evident that the pattern of the 15 correlations obtained by the simultaneous two-waveform fitting (AB in Fig. 6) is distinct from those of the single waveform fittings of the control AP (A)  or the I Kr blocked AP (B).Thus, it is clear that the PO method based on two-waveform fitting is provided with more information if compared to the single waveform fitting.The mechanisms underlying the variation of the correlation pattern may be complicated.Only, the disappearance of the correlation between I hasf and any others in panel B was attributed to I ha = 0.This is because the MDP is depolarized by the I Kr block to a voltage range, where I ha is totally deactivated.

Physiological ionic mechanisms underlying the optimized AP waveforms
So far, the application of the simultaneous two-waveform fitting of the PS method was examined at the methodological point if the ionic current profile expressed on the surface membrane was correctly determined by using the error function MSE of the PO method independently of the initial parameter set.It should also be asked if the profile of current expression levels obtained by the PO method well matches the biophysical function of the cell, which is initiated by the membrane excitation relying on the physiological characteristics of individual currents.
As demonstrated in Fig. 7A, the V m (black trace) closely followed the time course of V L (Eqs.15 and 16); V L > V m during the depolarizing phase and vice versa (V L < V m ) during the repolarizing phase.The V L crossed V m at dV m /dt = 0 at OS and MDP.Since the time course of V L was quite close to that of V m (Fig. 7A), here the contribution of each current to V m is described based on the corresponding changes in elements of both V L and dV L /dt in Fig. 8.
Figure 7B,C shows the time course of individual membrane currents underlying the AP waveform (V m ).In the control waveform in column 1 of Fig. 7 which corresponds to c in Fig. 3C, the I Kr (blue) rapidly showed a large peak at about − 50 mV and then slowly decayed during SDD.The time course of I NaK (pink) reflected its V m -dependency and the amplitude was relatively large if compared with others.I K1 (steel blue) was nearly zero during the plateau but showed a significant amplitude during SDD.The I Kur (purple) was activated only during AP and provided the largest current during AP.As to the inward currents, I bNSC (light brown) was the largest current during SDD and its amplitude was proportional to V m because of the absence of time-dependent gating kinetics in I bNSC .The I NCX current (red) transiently peaked in the outward direction at the foot of the AP upstroke, and the time course of the inward component reflected the Ca 2+ transient; the inward peak appeared ~ 50 ms after the rising phase of AP, and thereafter it decayed.During the late half of SDD, the amplitude of both I Na (purple) and I CaL (black) increased exponentially due to the subthreshold V m -dependent activation.I CaL was the major inward current during the plateau phase.In columns 2 and 3 of Fig. 7, I Kr decreased gradually from control in column 1.Although I xsf s of all current species other than I Kr were the same in three conditions in columns 1-3, the individual amplitude of each current in panels B and C varied according to their own kinetics when V m changed as the result of I Kr block.
The merit of the V L analysis is its quantitative assessment of the contribution of individual membrane currents in mV, or mV/ms unit (in the case of V L , or dV L /dt, respectively) to the AP waveform (Fig. 7).Indeed, the V L (gray) shown in Fig. 8A as well as in Fig. 7A (red) consists of the arithmetic sum of individual V L elements (Eq.15).Major V L elements were I NCX , I Kr , I K1 , I NaK and I bNSC during SDD, and I CaL , I NCX , and I NaK during AP plateau.It should be noted that the time courses of individual V L elements were all unique and quite different from those of the membrane currents, which provided a vital clue in exploring the direction of the PO by the PS method.
The time courses of individual elements of V L were determined by the changes in the channel open probability as well as by the difference in (V m − E x ) in the case of the channel current.The time course of I NaK and I NCX elements are dependent not only on V m but also largely determined by the [Na + ] i and [Ca 2+ ] i .The Ca 2+ -induced Ca 2+ -release (CICR) is initiated by the Ca 2+ influx through the LCC in the dyadic junction space.Then, the Ca 2+ -release channel on the SR membrane is activated by the massive increase in Ca 2+ in the dyadic junction.The time course of this CICR is clearly visualized in the dV L /dt records (Fig. 8B,C).Namely, the onset of I CaL activation is indicated by the rapid transient deflection of I CaL element of dV L /dt which coincides with the rise of I CaL element of V L .Then the transient upward deflection of the I NCX element of dV L /dt indicates that the driving force (reversal potential of I NCX − V m ) is reversed during the overshoot of AP.Then the relatively slow downward deflection of the I NCX element is driven by the increase of [Ca 2+ ] i within the dyadic space.This sequence of events was evident in the control (column 1 of Fig. 8) as well as in the two cases of I Kr -block (columns 2 and 3 of Fig. 8).The duration of SDD was shortened by both the delayed repolarization and the accelerated subthreshold activation of both I Na and I CaL because of the positive shift of MDP, both of which were caused by the I Kr -block (columns 2 and 3 of Fig. 8).With a more extensive I Kr block (column 3 of Fig. 8), the SDD virtually disappeared because of the overlap between the delayed repolarization on the subthreshold activation of I CaL and I Na .
These findings are totally consistent with the physiological mechanisms supposed so far through the experimental findings 1-4,10,40 .

Brief summary of results
(1) The parameter optimization method has been developed to determine the ionic current composition underlying the cardiac AP waveform.However, the method has not been established to be used in a wide range of electrophysiological research.It has been suggested that information included in a single AP waveform is still not enough for the PO method.In the present study, the AP record under the selective I Kr block was used in addition to the control AP waveform.(2) The MM test of the simultaneous fitting of two AP waveforms was applied to the AP waveform modified by blocking one of the currents, I Kr , I CaL , I Na , I Ks , I Kur , or I bNSC (Figs. 1 and 2, MM mode).In all cases of ion current block, the PS method successfully identified the blocked current, and the MSE was decreased fairly below log(MSE) < − 2. (3) In the EM test, the experimental record of the spontaneous AP frequently showed intermittent slow V m fluctuations, which interfered with the gradient-based PO method.However, by using the multi-run test, the distribution of the I xsf converged to a single point, even though the limit level of MSE decrease was largely enlarged in the EM test of the PS method (Fig. 4).The new method well resolved the I Kr -selective block current.(4) The correlations between two I xsf s were determined in 15 pairs of I xsf s.The pattern of 15 correlations was quite different among the three results obtained in the two-waveform fitting, and the single waveform fitting to the control and the I Kr block recordings (Fig. 6).This finding supports the view that the simultaneous fitting method is provided with more information than the single AP waveform fitting, and resolves the unique current pattern underlying the membrane excitation in the cell-specific computer model.(5) The current flows of the major 7 currents in the optimized cell-specific model were analyzed by the V L analysis method (Fig. 8).All these quantitative findings on the contributions of all currents were totally The importance of providing some information is required in addition to the single AP waveform when the PO method was applied to the cardiac AP.This point of improving the PO analysis has been extensively discussed in the literatures (see Kernik et al. 41 , Paci et al. 42 , Paci et al. 43 ).

The simultaneous fitting of two AP waveforms in the PS method.
The error function MSE is calculated from the difference between the target and model output AP waveforms to drive the PO.The waveform is determined by both the limiting conductance G x and time-dependent changes in pO as described by Eq. ( 3).The pO is a critical factor in shaping the precise AP waveform according to the detailed V m -dependent rate constant for the channel gating (Eqs.4 and 5).Particularly, the V m -dependency is quite different between different channels, so individual channel species can play distinct roles in different https://doi.org/10.1038/s41598-024-63413-0www.nature.com/scientificreports/

Figure 1 .
Figure 1.Convergence of I xsf (BP) as the MSE magnitude was decreased by the simultaneous two-waveform fitting in the MM mode multi-run test.The BPs were plotted in the coordinate of log(MSE) (ordinate) vs I xsf (abscissa) during the PO process by the PS method.In each line of graphs, I Krsf , I Nasf , I CaLsf , I Kssf , and I bNSCsf log(MSE) were plotted.In the case of the blocked current, the I xsf was normalized in reference to the blocked level.

Figure 2 .
Figure 2. The recovery of the AP configuration from the artificial modification of AP by the PO method.The red trace was the target AP, which was generated by the intact baseline model, while the black trace was the model output generated at the end of the PO method as the result.The violet horizontal lines indicate three ranges of AP to weight (6, 1, and 0.3 for each step, independent from the ordinate scale) in calculating the weighted sum of MSE.The fitting to the control trace in (A) was obtained by the simultaneous two-APwaveform fitting using the 60% I Kr block.Since the fittings of the control AP in cases of any other I x block were quite similar to that illustrated in (A), the respective control record for each I x block was not demonstrated individually.Note that two different percentages of blockage (30% and 60%) were specifically tested for I Kr .

Figure 3 .
Figure 3. Continuous recording of the spontaneous AP generation in hiPSC-CM during the I Kr block by perfusing E-4031.(A-C) Shows three examples of chart recordings of the AP waveform.The red traces indicate the single AP waveforms used for the two-waveform fitting of the PO method, which will be described in the next chapter.(D) Shows time courses of I xsf of I Kr (light blue), I CaL (orange), I Na (steel blue), I Kur (yellow) and I bNSC (gray) obtained by the conventional single waveform fitting of the PS method applied to TC12.See the text for a detailed explanation.

Figure 4 .
Figure 4. Results of the simultaneous fitting of two AP waveforms in the EM mode of the PO method.The data in three hiPSC-CMs (TC11, TC12, and TC13) were examined for the pair of control and the I Kr -blocked AP waveforms; the pair of c and b1, and the pair of c and b2 indicated in Fig. 3.In all panels the control APs (red) were superimposed by the model output APs (black).The results of the multi-run test are shown in the right half of figure.

Figure 5 .
Figure 5.Comparison of the I xsf s other than I Krsf examined in the EM test.The results of the two-waveform fitting using the leftmost and middle red traces shown in Fig. 3 are shown in the blue bars, and the results by using the leftmost and rightmost red traces shown in Fig. 3 are shown in the green bars.Note that the blue bars were normalized to 1.0 to compare with the green bars.See the text for the explanation.

Figure 6 .
Figure 6.The patterns of the correlations between different pairs of I xsf s were compared among three different POs.The correlation of all 15 pairs of I xsfs using the top 20 results obtained in the multi-run PO method for TC12 by using the leftmost and middle red traces shown in Fig. 3C (AB), by using the leftmost red trace only (A), and by using the middle red trace only (B) are shown with the lines of different thickness.The thicker lines indicate R 2 higher than 0.8, 0.6 and 0.4.See text for more details.

Figure 7 .
Figure 7.The profile of the current flow underlying the AP waveform of the cell-specific model of TC12 optimized by the PS method.The upper row (A) indicates V m (black) and V L (red).A1 indicates the control condition (c in Fig. 3C), and A2, A3 indicate the record obtained at two different degrees of the I Kr block (b1 and b2 in Fig. 3C).The lower rows (B) and (C) indicate outward and inward I m , respectively.The individual I m are depicted in different colors.The time scales indicated at the bottom are applied to both V m and I m panels.The yellow colors are used to highlight the mechanisms mainly for both the foot and plateau phases of AP.The peak of inward I CaL is far beyond the graphic display.

Figure 8 .
Figure 8.The profile of V L and dV L /dt elements of the cell-specific model of TC12 revealed by the V L analysis.The upper row (A) indicates V L elements.A1 indicates the control condition (c in Fig. 3C), and A2, A3 indicate the record obtained at two different degrees of the I Kr block.The lower rows (B) and (C) indicate dV L /dt elements for outward and inward currents, respectively.The individual V L and dV L /dt elements are depicted in different colors.The time scales indicated at the bottom are applied to both V L and dV L /dt records.The yellow colors are used to highlight the mechanisms underlying both the foot and plateau phases of AP.See the text for more explanation.

Table 1 .
List of the final I xsf and the [Na + ] i and [Ca] t obtained by the PS method.See the text for a more detailed explanation.The unit of [Na + ] i is mM, and that of [Ca] t is attomole.[Ca] t stands for the sum of ionized and bound forms of Ca (Ca tot ) in jnc, iz, blk, and SR.Vol indicates the volume of each compartment.These concentrations were measured at the time when Vm crosses − 20 mV during the repolarization phase.Ca tot = Ca tot,jnc • Vol jnc + Ca tot,iz • Vol iz + Ca tot,blk • Vol blk + Ca SRup • Vol SRup + Ca tot,SRrl • Vol SRrl .